<?php

/**
 * @package JAMA
 *
 *    For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
 *    orthogonal matrix Q and an n-by-n upper triangular matrix R so that
 *    A = Q*R.
 *
 *    The QR decompostion always exists, even if the matrix does not have
 *    full rank, so the constructor will never fail.  The primary use of the
 *    QR decomposition is in the least squares solution of nonsquare systems
 *    of simultaneous linear equations.  This will fail if isFullRank()
 *    returns false.
 *
 * @author  Paul Meagher
 * @license PHP v3.0
 * @version 1.1
 */
class PHPExcel_Shared_JAMA_QRDecomposition {

	const MatrixRankException = "Can only perform operation on full-rank matrix.";

	/**
	 *    Array for internal storage of decomposition.
	 *
	 * @var array
	 */
	private $QR = array();

	/**
	 *    Row dimension.
	 *
	 * @var integer
	 */
	private $m;

	/**
	 *    Column dimension.
	 *
	 * @var integer
	 */
	private $n;

	/**
	 *    Array for internal storage of diagonal of R.
	 *
	 * @var  array
	 */
	private $Rdiag = array();


	/**
	 *    QR Decomposition computed by Householder reflections.
	 *
	 * @param matrix $A Rectangular matrix
	 *
	 * @return Structure to access R and the Householder vectors and compute Q.
	 */
	public function __construct($A) {
		if ($A instanceof PHPExcel_Shared_JAMA_Matrix)
		{
			// Initialize.
			$this->QR = $A->getArrayCopy();
			$this->m = $A->getRowDimension();
			$this->n = $A->getColumnDimension();
			// Main loop.
			for ($k = 0; $k < $this->n; ++$k)
			{
				// Compute 2-norm of k-th column without under/overflow.
				$nrm = 0.0;
				for ($i = $k; $i < $this->m; ++$i)
				{
					$nrm = hypo($nrm, $this->QR[$i][$k]);
				}
				if ($nrm != 0.0)
				{
					// Form k-th Householder vector.
					if ($this->QR[$k][$k] < 0)
					{
						$nrm = -$nrm;
					}
					for ($i = $k; $i < $this->m; ++$i)
					{
						$this->QR[$i][$k] /= $nrm;
					}
					$this->QR[$k][$k] += 1.0;
					// Apply transformation to remaining columns.
					for ($j = $k + 1; $j < $this->n; ++$j)
					{
						$s = 0.0;
						for ($i = $k; $i < $this->m; ++$i)
						{
							$s += $this->QR[$i][$k] * $this->QR[$i][$j];
						}
						$s = -$s / $this->QR[$k][$k];
						for ($i = $k; $i < $this->m; ++$i)
						{
							$this->QR[$i][$j] += $s * $this->QR[$i][$k];
						}
					}
				}
				$this->Rdiag[$k] = -$nrm;
			}
		}
		else
		{
			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
		}
	} //	function __construct()


	/**
	 *    Is the matrix full rank?
	 *
	 * @return boolean true if R, and hence A, has full rank, else false.
	 */
	public function isFullRank() {
		for ($j = 0; $j < $this->n; ++$j)
		{
			if ($this->Rdiag[$j] == 0)
			{
				return false;
			}
		}

		return true;
	} //	function isFullRank()


	/**
	 *    Return the Householder vectors
	 *
	 * @return Matrix Lower trapezoidal matrix whose columns define the reflections
	 */
	public function getH() {
		for ($i = 0; $i < $this->m; ++$i)
		{
			for ($j = 0; $j < $this->n; ++$j)
			{
				if ($i >= $j)
				{
					$H[$i][$j] = $this->QR[$i][$j];
				}
				else
				{
					$H[$i][$j] = 0.0;
				}
			}
		}

		return new PHPExcel_Shared_JAMA_Matrix($H);
	} //	function getH()


	/**
	 *    Return the upper triangular factor
	 *
	 * @return Matrix upper triangular factor
	 */
	public function getR() {
		for ($i = 0; $i < $this->n; ++$i)
		{
			for ($j = 0; $j < $this->n; ++$j)
			{
				if ($i < $j)
				{
					$R[$i][$j] = $this->QR[$i][$j];
				}
				elseif ($i == $j)
				{
					$R[$i][$j] = $this->Rdiag[$i];
				}
				else
				{
					$R[$i][$j] = 0.0;
				}
			}
		}

		return new PHPExcel_Shared_JAMA_Matrix($R);
	} //	function getR()


	/**
	 *    Generate and return the (economy-sized) orthogonal factor
	 *
	 * @return Matrix orthogonal factor
	 */
	public function getQ() {
		for ($k = $this->n - 1; $k >= 0; --$k)
		{
			for ($i = 0; $i < $this->m; ++$i)
			{
				$Q[$i][$k] = 0.0;
			}
			$Q[$k][$k] = 1.0;
			for ($j = $k; $j < $this->n; ++$j)
			{
				if ($this->QR[$k][$k] != 0)
				{
					$s = 0.0;
					for ($i = $k; $i < $this->m; ++$i)
					{
						$s += $this->QR[$i][$k] * $Q[$i][$j];
					}
					$s = -$s / $this->QR[$k][$k];
					for ($i = $k; $i < $this->m; ++$i)
					{
						$Q[$i][$j] += $s * $this->QR[$i][$k];
					}
				}
			}
		}

		/*
		for($i = 0; $i < count($Q); ++$i) {
			for($j = 0; $j < count($Q); ++$j) {
				if(! isset($Q[$i][$j]) ) {
					$Q[$i][$j] = 0;
				}
			}
		}
		*/

		return new PHPExcel_Shared_JAMA_Matrix($Q);
	} //	function getQ()


	/**
	 *    Least squares solution of A*X = B
	 *
	 * @param Matrix $B A Matrix with as many rows as A and any number of columns.
	 *
	 * @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
	 */
	public function solve($B) {
		if ($B->getRowDimension() == $this->m)
		{
			if ($this->isFullRank())
			{
				// Copy right hand side
				$nx = $B->getColumnDimension();
				$X = $B->getArrayCopy();
				// Compute Y = transpose(Q)*B
				for ($k = 0; $k < $this->n; ++$k)
				{
					for ($j = 0; $j < $nx; ++$j)
					{
						$s = 0.0;
						for ($i = $k; $i < $this->m; ++$i)
						{
							$s += $this->QR[$i][$k] * $X[$i][$j];
						}
						$s = -$s / $this->QR[$k][$k];
						for ($i = $k; $i < $this->m; ++$i)
						{
							$X[$i][$j] += $s * $this->QR[$i][$k];
						}
					}
				}
				// Solve R*X = Y;
				for ($k = $this->n - 1; $k >= 0; --$k)
				{
					for ($j = 0; $j < $nx; ++$j)
					{
						$X[$k][$j] /= $this->Rdiag[$k];
					}
					for ($i = 0; $i < $k; ++$i)
					{
						for ($j = 0; $j < $nx; ++$j)
						{
							$X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
						}
					}
				}
				$X = new PHPExcel_Shared_JAMA_Matrix($X);

				return ($X->getMatrix(0, $this->n - 1, 0, $nx));
			}
			else
			{
				throw new PHPExcel_Calculation_Exception(self::MatrixRankException);
			}
		}
		else
		{
			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
		}
	} //	function solve()

} //	PHPExcel_Shared_JAMA_class QRDecomposition
